System and method for volumetric tumor segmentation using joint space-intensity likelihood ratio test

ABSTRACT

A method for segmenting a digitized image includes providing a digitized volumetric image comprising a plurality of intensities corresponding to a domain of points in an N-dimensional space, identifying a target structure in said image, forming a window about said target structure whose size is a function of the target scale, and performing a joint space-intensity-likelihood ratio test at each point within said window to determine whether each said point is within said target structure.

CROSS REFERENCE TO RELATED UNITED STATES APPLICATIONS

This application claims priority from “Volumetric Tumor Segmentationusing Space-Intensity Joint Likelihood Ratio Test”, U.S. ProvisionalApplication No. 60/608,499 of Okada, et al., filed Sep. 9, 2004, andfrom “Blob Segmentation using Joint Space-Intensity Likelihood RatioTest: Application to 3D Tumor Segmentation”, U.S. ProvisionalApplication No. 60/625,027 of Okada, et al., filed Nov. 4, 2004, thecontents of both of which are incorporated herein by reference.

TECHNICAL FIELD

This invention is directed to object segmentation in digitized medicalimages.

DISCUSSION OF THE RELATED ART

The diagnostically superior information available from data acquiredfrom current imaging systems enables the detection of potential problemsat earlier and more treatable stages. Given the vast quantity ofdetailed data acquirable from imaging systems, various algorithms mustbe developed to efficiently and accurately process image data. With theaid of computers, advances in image processing are generally performedon digital or digitized images.

Digital images are created from an array of numerical valuesrepresenting a property (such as a grey scale value or magnetic fieldstrength) associable with an anatomical location points referenced by aparticular array location. The set of anatomical location pointscomprises the domain of the image. In 2-D digital images, or slicesections, the discrete array locations are termed pixels.Three-dimensional digital images can be constructed from stacked slicesections through various construction techniques known in the art. The3-D images are made up of discrete volume elements, also referred to asvoxels, composed of pixels from the 2-D images. The pixel or voxelproperties can be processed to ascertain various properties about theanatomy of a patient associated with such pixels or voxels.

The process of classifying, identifying, and characterizing imagestructures is known as segmentation. Once anatomical regions andstructures are identified by analyzing pixels and/or voxels, subsequentprocessing and analysis exploiting regional characteristics and featurescan be applied to relevant areas, thus improving both accuracy andefficiency of the imaging system. The wide variety of object appearancecharacteristics and boundary geometry makes image segmentation a verydifficult task. In past decades, a number of promising general-purposeapproaches, such as classification/labeling/clustering andcurve-evolution, have been proposed to solve this problem. In practice,however, structural assumptions of the target objects are oftenavailable beforehand thus can be exploited as a prior. The successfulincorporation of such prior information plays a key role for realizingefficient and accurate segmentation solutions in general.

The development of medical data segmentation solutions as applied tocomputer aided diagnosis applications emphasizes the overall systemperformance, including user-interaction factors. In such context,semi-automatic solutions, requiring minimal user interactions, can bepreferred to fully automated solutions for achieving better overallperformance. For this reason, a one-click figure-ground segmentationapproach is preferred, where a user can provide a data point whichroughly indicates a target/figure blob to be segmented out of arbitrarybackground. A successful solution depends on (1) robustness againstvariation of the user-given initialization and the different scansettings to relieve the user's labor, (2) run-time efficiency, even withthe high-dimensional data, to enhance the user-interactivity, and (3)high accuracy so that the user-interaction results in better performancethan a fully automated solution.

SUMMARY OF THE INVENTION

Exemplary embodiments of the invention as described herein generallyinclude methods and systems for semi-automatic figure-groundsegmentation solution for blob-like objects in multi-dimensional images.The blob-like structure include various objects of interest that arehard to segment in many application domains, such as tumor lesions in 3Dmedical data. The embodiment of the present invention are motivatedtowards computer-aided diagnosis medical applications, justifying asemi-automatic figure-ground approach. An efficient segmentation isrealized by combining anisotropic Gaussian model fitting and alikelihood ratio test (LRT)-based nonparametric segmentation in jointspace-intensity domain. The robustly fitted Gaussian is exploited toestimate the foreground and background likelihoods for both spatial andintensity variables. The LRT with the bootstrapped likelihoods is theoptimal Bayesian classification while automatically determining the LRTthreshold. A 3D implementation of one embodiment is applied to the lungnodule segmentation in CT data and validated with 1310 cases. A targetnodule is segmented in less than 3 seconds in average.

According to an aspect of the invention, there is provided a method forsegmenting a digitized image comprising the steps of providing adigitized volumetric image comprising a plurality of intensitiescorresponding to a domain of points in an N-dimensional space, providingan approximate location of a target structure in said image, estimatinga foreground spatial-intensity likelihood function about said targetstructure, estimating a background spatial-intensity likelihood functionabout said target structure, and using said foreground and backgroundspatial-intensity likelihood functions to segment said target structureby determining whether a point about said target structure is insidesaid target structure.

According to a further aspect of the invention, the method comprisesdetermining an estimated center and an estimated spread of said targetstructure by fitting an N-dimensional anisotropic Gaussian function to avolume of interest centered about the approximate location anddetermining the center and the anisotropic spread of said Gaussianfunction.

According to a further aspect of the invention, the foregroundspatial-intensity likelihood function can be factored into a product ofa foreground spatial likelihood function and a foreground intensitylikelihood function, and said background spatial-intensity likelihoodfunction can be factored into a product of a background spatiallikelihood function and a background intensity likelihood function.

According to a further aspect of the invention, the foreground spatiallikelihood function is proportional to said anisotropic Gaussianfunction, and said background intensity likelihood function is acomplement of said foreground spatial likelihood function.

According to a further aspect of the invention, the method comprisesimposing a window about said target structure, wherein said window isdefined as those points whose Mahalanobis distance from said mean ofsaid Gaussian is less than a predetermined constant value, wherein saidMahalanobis distance is computed using said spread of said Gaussian.

According to a further aspect of the invention, the constant value isdetermined by solving${{{{2{\pi\Sigma}}}^{{- 1}/2}{\int_{S{(c)}}{{\exp\left( {{- \frac{1}{2}}\left( {x - u} \right)^{t}{\sum\limits^{- 1}\left( {x - u} \right)}} \right)}{\mathbb{d}x}}}} = \frac{{S(c)}}{2{{2{\pi\Sigma}}}^{1/2}}},$wherein Σ is said spread, c is said constant value, S(c) is said window,x is a point in said window, and u is the center of said targetstructure, and S(c) = ∫_(S(c))𝕕x.

According to a further aspect of the invention, the foreground intensitylikelihood function is proportional to a foreground intensity differencefunction weighted by said foreground spatial likelihood function sampledwithin said window, and said background intensity likelihood function isproportional to a background intensity difference function weighted bysaid background spatial likelihood function sampled within said window.According to a further aspect of the invention, the proportionalityconstant is equal to one half the norm of the window.

According to a further aspect of the invention, the foreground andbackground intensity difference functions comprise Dirac deltafunctions.

According to a further aspect of the invention, the foreground andbackground intensity difference functions comprise Parzen functions.

According to a further aspect of the invention, the step of determiningwhether a point about said target structure is inside said targetstructure is repeated for every point neighboring said target structureto determine which points comprise said target structure.

According to a further aspect of the invention, determining whether apoint is inside said target structure comprises comparing a ratio ofsaid foreground and background spatial-intensity likelihood functionscalculated at said point to a preset threshold, wherein said point isclassified as inside said target structure if said ratio is greater thansaid threshold.

According to a further aspect of the invention, determining whether apoint is inside said target structure comprises comparing saidforeground spatial-intensity likelihood function ƒ(x, α|in) to saidbackground spatial-intensity likelihood function ƒ(x, α|out), whereinsaid point x is classified as inside said target structure if ƒ(x,α|in)>ƒ(x, α|out).

According to a further aspect of the invention, determining whether apoint is inside said target structure comprises comparing a function Fof a ratio of said foreground likelihood function to said backgroundlikelihood function to F(1) at the point location, wherein the functionF is a member of a family of functions F:R→R that are monotonically andstrictly increasing, wherein said point is classified as inside saidtarget structure if the function of said ratio is greater than F(1).

According to a further aspect of the invention, determining whether apoint is inside said target structure comprises comparing a function Fof said foreground likelihood function ƒ(x, α|in) to a function F ofsaid background likelihood function ƒ(x, α|out) at the point location x,wherein the function F is a member of a family of functions F:R→R thatare monotonically and strictly increasing, and wherein the point x isclassified as inside said target structure if F(ƒ(x, α|in))≧F(ƒ(x,α|out).

According to another aspect of the invention, there is provided aprogram storage device readable by a computer, tangibly embodying aprogram of instructions executable by the computer to perform the methodsteps for segmenting a digitized image.

BRIEF DESCRIPTION OF THE DRAWINGS

FIGS. 1(a)-(f) illustrate the likelihood estimation processes for a 1Dexample, according to an embodiment of the invention.

FIG. 2 depicts a flow chart of a joint space-intensity likelihood ratiotest based segmentation method, according to an embodiment of theinvention.

FIGS. 3(a)-(d) illustrate examples of segmentation results for fourtumor cases, according to an embodiment of the invention.

FIGS. 4(a)-(d) shows the intensity likelihood models estimated for thefour cases in FIGS. 3(a)-(d), according to an embodiment of theinvention.

FIG. 5 illustrates examples of 2D views of 3D segmentation results forfive tumor cases, according to an embodiment of the invention.

FIG. 6 is a block diagram of an exemplary computer system forimplementing a joint space-intensity likelihood ratio test basedsegmentation method, according to an embodiment of the invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

Exemplary embodiments of the invention as described herein generallyinclude systems and methods for an efficient segmentation solution for aclass of blob-like structures captured in multi-dimensional medicalimages. Although an exemplary embodiment of this invention is discussedin the context of segmenting a CT lung nodule, it is to be understoodthat the object segmentation and shape characterization methodspresented herein have application to other multi-dimensional imagingmodalities.

As used herein, the term “image” refers to multi-dimensional datacomposed of discrete image elements (e.g., pixels for 2-D images andvoxels for 3-D images). The image may be, for example, a medical imageof a subject collected by computer tomography, magnetic resonanceimaging, ultrasound, or any other medical imaging system known to one ofskill in the art. The image may also be provided from non-medicalcontexts, such as, for example, remote sensing systems, electronmicroscopy, etc. Although an image can be thought of as a function fromR³ to R, the methods of the inventions are not limited to such images,and can be applied to images of any dimension, e.g. a 2-D picture or a3-D volume. For a 2- or 3-dimensional image, the domain of the image istypically a 2- or 3-dimensional rectangular array, wherein each pixel orvoxel can be addressed with reference to a set of 2 or 3 mutuallyorthogonal axes. The terms “digital” and “digitized” as used herein willrefer to images or volumes, as appropriate, in a digital or digitizedformat acquired via a digital acquisition system or via conversion froman analog image.

A blob-like structure can be defined as a roughly convex local intensitydistribution whose iso-level contours are approximately ellipsoidal, butwith some irregularities that do not destroy the ellipsoidal topology.The intensity distribution itself can be multi-modal but can be assumedto be uni-modal under Gaussian blurring within an appropriateupper-bound of the smoothing bandwidth. Such a class of data structuresrepresents various objects of interest, such as tumors and polyps, whichare hard to segment in many medical imaging application domains, such asCT lung and PET hot spot segmentation.

According to an embodiment of the invention, a semi-automatic(one-click) blob segmentation method includes two steps. An example of ablob is a tumor. A first step is a pre-processing step with aanisotropic Gaussian fitting. Given an initial marker x indicating anapproximate location of the target structure, such as a tumor, theGaussian fitting provides an estimated target center u and ananisotropic spread matrix Σ in the form of Gaussian function:${\Phi\left( {x,u,\Sigma} \right)} = {\frac{\exp\left( {{- \frac{1}{2}}\left( {x - u} \right)^{t}{\sum\limits^{- 1}\left( {x - u} \right)}} \right)}{{{2{\pi\Sigma}}}^{1/2}}.}$Note that the notation (. . . )^(t) indicated the transpose of a vector(or matrix). The volume of interest (VOI) Ω can be defined by the extentof the data analysis given by a fixed-size N-D window centered at x. Thedata to be analyzed is expressed by I(x)εR⁺ where xεΩ⊂R^(N) is anN-dimensional coordinate indicating a data (pixel/voxel) location. Theresulting multi-scale Gaussian model fitting solution is robust against(1) the influence from non-target neighboring structures, (2) misfit ofthe data, and (3) variations in the initialization point x. Theanisotropic Gaussian fitting procedure is described in the inventors'copending patent application, “Method for Robust Scale-Space Analysis of3D Local Structures in Medical Images”, U.S. patent application Ser. No.10/892,646, filed Jul. 17, 2004, the contents of which are incorporatedherein by reference.

A second step of the segmentation method according to an embodiment ofthe invention involves using a likelihood test to separate a figure fromthe background. At each data point xεΩ there is an intensity valueα=I(x). Treating both x and α as independent random variables, a jointlikelihood function of (x, α) can be estimated for a foreground, ƒ(x,α|in), where in represents the interior of or part of a target tumor,and for a background, ƒ(x, α|out), where out represents the outside ofthe tumor. The space-intensity joint likelihoods can be factorized asƒ(x, α|in)=ƒ(x|in)ƒ(α|in)ƒ(x, α|out)=ƒ(x, out)ƒ(α|out)where ƒ(x|in) and ƒ(α|in), ƒ(x|out) and ƒ(α|out)) denote a marginalforeground (background) spatial and intensity likelihood functions,respectively. Although the two variables x and α are not independent ingeneral, experimental results have shown indicate that their dependenceis weak, resulting in good segmentation results. A space-intensity jointlikelihood ratio r(x) is then defined by${{r(x)} \equiv \frac{f\left( {x,{\alpha ❘{in}}} \right)}{f\left( {x,{\alpha ❘{out}}} \right)}} = {\frac{{f\left( {x❘{in}} \right)}{f\left( {\alpha ❘{in}} \right)}}{{f\left( {x❘{out}} \right)}{f\left( {\alpha ❘{out}} \right)}}.}$Each voxel data point within the VOI can be segmented by performing thelikelihood ratio test: xεin if r(x)≧th, otherwise xεout, where th is athreshold which depends on the normalization factors of the foregroundand background likelihoods. Modeling the likelihoods within a specificsupport region assures the Bayesian optimality for th=1. It is to beunderstood, however, that other threshold values can be used, and thethreshold can vary for different sets of points.

It is to be further understood that the likelihood ratio is an exemplaryfunction of the foreground and background joint space-intensitylikelihood functions, and other tests involving these likelihoodfunctions can be used and be within the scope of an embodiment of theinvention. In one exemplary embodiment, likelihood-based segmentation isperformed by comparison of the foreground and background jointspace-intensity likelihood functions. Given the positive-valuedforeground likelihood function ƒ(x, α|in) and the background likelihoodfunction ƒ(x, α|out) at a point location x, the point x is classified asa member of the foreground if ƒ(x, α|in)>ƒ(x, α|out), otherwise, it isclassified as the background. This segmentation scheme is equivalent tothe likelihood ratio test based segmentation method. It is preferred tothe likelihood ratio test when the background likelihood ƒ(x, α|out) canbe zero-valued for some locations x, where the foreground/backgroundration would not be computable.

In another exemplary embodiment, given the positive-valued foregroundlikelihood function ƒ(x, α|in) and the background likelihood functionƒ(x, α|out) at a point location x, the point x is classified as a memberof the foreground if some function F(ƒ(x, α|in)/ƒ(x, α|out)) is largerthan or equivalent to F(1) and classified as the background otherwise.The function F is a member of a family of functions F:R→R that aremonotonically and strictly increasing, i.e., order-preserving.

In another embodiment of the invention, given the positive-valuedforeground likelihood function ƒ(x, α|in) and the background likelihoodfunction ƒ(x, α|out) at the point location x, the point x is classifiedas a member of the foreground if F(ƒ(x, α|in))≧F(ƒ(x, α|out) andclassified as the background otherwise. Again, the function F is amember of a family of functions F:R→R that are monotonically andstrictly increasing. Examples of functions that are monotonically andstrictly increasing include logarithmic functions, polynomial functions,and exponential functions.

However, it is to be understood that not all functions of the likelihoodfunctions will yield the desired results. Although tests such as ƒ(x,α|in)>ƒ(x, α|out), log(ƒ(x, α|in))>log(ƒ(x, α|out)), sqrt(ƒ(x,α|in)/ƒ(x, α|out))>1, etc., will yield consistent results, other tests,such as (ƒ(x, α|in))²>ƒ(x, α|out) and log(ƒ(x, α|in))/log(ƒ(x,α|out))>th, will not yield results consistent with the likelihood ratiotest method ƒ(x, α|in)/ƒ(x, α|out)>1.

One can realize the segmentation outlined above by defining fourlikelihood functions for spatial and intensity factors inside andoutside the target structure: ƒ(x|in), ƒ(x|out),ƒ(α|in), and ƒ(α|out).

One can obtain the foreground and background spatial likelihoods byassuming that the N-D Gaussian model fitting solution function Φ(x, u,Σ) approximates a probability distribution of a location x being a tumorcenter or mean u. In many applications, such as tumor segmentation, thesurface geometry of the target structure is approximately convex, whichassures the mean is located inside of the structure. Thus, according toan embodiment of the invention, the Gaussian model fitting solutionfunction can be interpreted as the conditional probability distributionP(x|in) of x being part of the target tumor structure: P(x|in)=Φ(x, u,Σ). However, the conditional probability distribution for thebackground, P(x|out) is ill-defined because the background has aninfinite extent in the data space x. According to another embodiment ofthe invention, a support window S⊂Ω that confines observations of therandom variable x can be introduced so that the background has a finitenormalization. A pair of normalized conditional probabilitydistributions functions can be defined over the support window as{overscore (P)}(x|in)≡P(x|in)/∫_(S) P(x|in)dx{overscore (P)}(x|out)≡P(x|out)/∫_(S) P(x|out)dxwhere P(x|in) is known and P(x|out) is an unknown underlying backgrounddistribution. The total probability of a data point being in S isP_(x)={overscore (P)}(x|in)P_(in)+{overscore (P)}(x|out)P_(out)=1|S|,where Pin and Pout are prior probabilities of being inside and outsideof S, P_(in)+P_(out)=1. Thus, according to an embodiment of theinvention, the background spatial probability distribution can bedefined as${\overset{\_}{P}\left( {x❘{out}} \right)} = {\frac{\left( {1/{s}} \right) - {{\overset{\_}{P}\left( {x❘{in}} \right)} \times P_{i\quad n}}}{\left( {1 - P_{i\quad n}} \right)}.}$According to another embodiment of the invention, the inside and outsideprobabilities are unbiased and can be set equal so thatP_(in)=P_(out)=0.5:${\overset{\_}{P}\left( {x❘{out}} \right)} = {{\frac{2}{S} - {\overset{\_}{P}\left( {x❘{in}} \right)}} = {\frac{2}{S} - \frac{P\left( {x❘{in}} \right)}{\int_{S}{{P\left( {x❘{in}} \right)}{\mathbb{d}x}}}}}$According to another embodiment of the invention, the backgroundprobability distribution function over S can assume the value zero atthe mean location u where the Gaussian function Φ(x, u, Σ) modeling theforeground probability distribution function, takes its maximum. In thisembodiment, the normalization factor of {overscore (P)}(x|in) becomes${{\int_{S}{{P\left( {x❘{in}} \right)}{\mathbb{d}x}}} = {{{P\left( {u❘{in}} \right)}\frac{|S|}{2}} = \frac{S}{2{{2{\pi\Sigma}}}^{1/2}}}},$and the normalized foreground and background distributions can bedefined as${\overset{\_}{P}\left( {x❘{in}} \right)} = {\frac{2}{S}{{2{\pi\Sigma}}}^{1/2}{P\left( {x❘{in}} \right)}}$${\overset{\_}{P}\left( {x❘{out}} \right)} = {\frac{2}{S}\left( {1 - {{{2{\pi\Sigma}}}^{1/2}{P\left( {x❘{in}} \right)}}} \right)}$The foreground and background spatial likelihood functions can bedefined in terms as the conditional probability distribution functionsover S scaled by a fixed factor |S|/2 so that they depend only onP(x|in):${{f\left( {x❘{in}} \right)} \equiv {\frac{S}{2}{\overset{\_}{P}\left( {x❘{in}} \right)}}} = {{{2{\pi\Sigma}}}^{12}{P\left( {x❘{in}} \right)}}$${{f\left( {x❘{out}} \right)} \equiv {\frac{S}{2}{\overset{\_}{P}\left( {x❘{out}} \right)}}} = {1 - {{{2{\pi\Sigma}}}^{1/2}{{P\left( {x❘{in}} \right)}.}}}$Note that the background likelihood ƒ(x|out) is a complement of theforeground likelihood. At the mean location u, ƒ(u|in)=1 and ƒ(u|out)=0.At infinity, ƒ(±≡|in)=0 and ƒ(±≡|out)=1. In addition, since thelikelihood functions share a common scaling factor, the ratio of thelikelihood functions is equivalent to the ratio of the distributionfunctions.

The choice of the support window S can effect the segmentation solution.As previously described, the background can have an infinite spatialextent, and thus a background spatial likelihood function is not boundedand would have an infinite normalization factor. For this reason, asupport window S was introduced so that probability distributions can bedefined within such a window. However, the estimated backgroundlikelihood will be sensitive to the varying range of S since suchvariation of the support S would cause a large change to thenormalization factor.

According to an embodiment of the invention, the support window S can bea function of the target scale. For example, if a cup on a table is tobe segmented, it is sensible to model the background using specificinformation from the table, not of the house in which the table isplaced nor of the city in which the house is present. The Gaussianfunction fitted to the target structure by the pre-process can providesuch scale information in the form of a confidence ellipsoid ofN-dimensional equal-probability contour approximating the structureboundary. Utilizing this, the support window S can be parameterized as afunction of the ellipsoid:S(c)≡{x|(x−u)^(t) Σ ⁻¹(x−u)≦c}where the scalar c is the Mahalanobis distance of x from u withcovariance Σ. The constant c can be determined from S(c) = ∫_(S(c))𝕕xand the normalization of P(x|in):${\int_{S{(c)}}{{P\left( {x❘{in}} \right)}{\mathbb{d}x}}} = {\frac{{S(c)}}{2{{2{\pi\Sigma}}}^{1/2}} = {{{2{\pi\Sigma}}}^{{- 1}/2}{\int_{S{(c)}}{{\exp\left( {{- \frac{1}{2}}\left( {x - u} \right)^{t}{\sum\limits^{- 1}\left( {x - u} \right)}} \right)}{{\mathbb{d}x}.}}}}}$The solution S(c) depends on the dimensionality N of the data space x.For example, numerical solutions of the above equation for 1D, 2D and 3Dcases are: c₁≈6.1152, c₂≈3.1871, and c₃≈2.4931. Within this supportwindow, the probability mass of ƒ(x|in) and ƒ(x|out) over S areequivalent.

For the 3D segmentation, c₃=2.4931 amounts to an approximate 52%confidence interval of the chi-square distribution with three degrees offreedom. Empirically, previous studies for 3D tumor segmentationindicate that the equal-probability contour with c₃=1.6416, derived from35% confidence interval of the fitted Gaussian function, approximatesthe tumor boundary well. This suggests that S(c₃) derived above providesa data range that covers the complete foreground and includes only athin layer of background region around the target. This is anappropriate support window for modeling the background because thebackground model estimated over this support window will not be stronglyinfluenced by the non-target neighboring structures that may appearwithin the VOI.

One can obtain the foreground and background intensity likelihoods bydefining the conditional intensity probability distributions as afunction of intensity differences weighted by the correspondingnormalized spatial probability distributions and sampled over thesupport window S:${{\overset{\_}{P}\left( {\alpha ❘{in}} \right)} \equiv {\int_{S}{{\overset{\_}{P}\left( {x,{\alpha ❘{in}}} \right)}{\mathbb{d}x}}}} = {\int_{S}{{\overset{\_}{P}\left( {x❘{in}} \right)}\phi\quad\left( {{I(x)} - \alpha} \right){\mathbb{d}x}}}$${{\overset{\_}{P}\left( {\alpha ❘{out}} \right)} \equiv {\int_{S}{{\overset{\_}{P}\left( {x,{\alpha ❘{out}}} \right)}{\mathbb{d}x}}}} = {\int_{S}{{\overset{\_}{P}\left( {x❘{out}} \right)}\phi\quad\left( {{I(x)} - \alpha} \right){\mathbb{d}x}}}$where {overscore (P)}(α|x,{in/out}) is modeled by φ(I(x)−α). Thefunction φ should be localized and have finite normalization. There areseveral possibilities for this function. In one embodiment of theinvention, to assure unit-normalization over the support window S, thefunction φ can be set to the discrete Dirac delta function. According toanother embodiment of the invention, for estimating a continuousprobability distribution function from a small number of samples, aParzen function with a uniform step kernel can be used as φ whilemaintaining the unit-normalization. Replacing the spatial conditionaldistributions by the likelihood functions yields${\overset{\_}{P}\left( {\alpha ❘{in}} \right)} = {\frac{2}{S}{\int_{S}{{f\left( {x❘{in}} \right)}{\phi\left( {{I(x)} - \alpha} \right)}{\mathbb{d}x}}}}$${\overset{\_}{P}\left( {\alpha ❘{out}} \right)} = {\frac{2}{S}{\int_{S}{{f\left( {x❘{out}} \right)}{\phi\left( {{I(x)} - \alpha} \right)}{\mathbb{d}x}}}}$Similar to the spatial likelihood functions, the intensity likelihoodfunctions can be defined as scaled conditional distribution functionswith a fixed factor |S|/2 sampled over the support window S:${{f\left( {\alpha ❘{in}} \right)} \equiv {\frac{S}{2}{\overset{\_}{P}\left( {\alpha ❘{in}} \right)}}} = {\int_{S}{{f\left( {x❘{in}} \right)}{\phi\left( {{I(x)} - \alpha} \right)}{\mathbb{d}x}}}$${{f\left( {\alpha ❘{out}} \right)} \equiv {\frac{S}{2}{\overset{\_}{P}\left( {\alpha ❘{out}} \right)}}} = {\int_{S}{{f\left( {x❘{out}} \right)}{\phi\left( {{I(x)} - \alpha} \right)}{\mathbb{d}x}}}$The techniques for likelihood estimation according to embodiments of thepresent invention do not require iterative model updates since theGaussian fitting step provides a robust and accurate targetcharacterization, captured in ƒ(x|in) and ƒ(x|out).

FIGS. 1(a)-(f) illustrate the likelihood estimation processes for a 1Dexample, according to an embodiment of the invention. FIG. 1(a) depicts1D noisy data with a Gaussian fitted by the pre-process. FIG. 1(b)depicts foreground (solid) and background (dash) spatial likelihoods andthe support window (dot-dash), derived from the Gaussian. Given thefitted Gaussian shown by a dash-curve in FIG. 1(a), the foreground(solid) and background (dash) likelihoods are analytically determined asshown in FIG. 1(b). Because both fore- and back-ground likelihoods sharethe same scaling factor, the ratio of the likelihoods and that of theprobability distribution functions become equivalent.

FIGS. 1(c)-(f) illustrate the intensity likelihood estimation processes.FIG. 1(c) illustrates the data, indicating with the solid lines a pairof pixel location and intensity values (x_(i), α_(i)). FIG. 1(d)illustrates the spatial likelihoods, showing the foreground (solid) andbackground (dashed) likelihoods at x_(i). FIG. 1(e) depicts theforeground intensity likelihood, showing contributions (dashed lines)from the data point (x_(i), α_(i)). FIG. 1(f) depicts the backgroundintensity likelihood, showing contribution (dashed lines) from (x_(i),α_(i)). Using all data within the support window (x_(i)εS, α_(i)), theforeground (FIG. 1(e)) and background (FIG. 1(f)) intensity likelihoodsare estimated by accumulating φ-smoothed counts for each intensity valueα_(i) weighted by the corresponding spatial likelihoods ƒ(x_(i)|in) andƒ(x_(i)|out) shown in FIG. 1(d).

With the spatial and intensity likelihood functions derived above, thejoint likelihood ratio r(x) can be expressed as:${r(x)} = \frac{{{2{\pi\Sigma}}}^{1/2}{\Phi\left( {x,u,\Sigma} \right)}{\int_{S}{{{2{\pi\Sigma}}}^{1/2}{\Phi\left( {x,u,\Sigma} \right)}{\phi\left( {{I(x)} - \alpha} \right)}{\mathbb{d}x}}}}{\left( {1 - {{{2{\pi\Sigma}}}^{1/2}{\Phi\left( {x,u,\Sigma} \right)}}} \right){\int_{S}{\left( {1 - {{{2{\pi\Sigma}}}^{1/2}{\Phi\left( {x,u,\Sigma} \right)}}} \right){\phi\left( {{I(x)} - \alpha} \right)}{\mathbb{d}x}}}}$This shows that the likelihood ratio at x with intensity value a dependsonly on Φ(x, u, Σ) and I(xεS). The formal derivations presented aboveassure that the ratios of the foreground and background likelihoods areequivalent to the ratios of the posterior probability distributionfunctions normalized over the support window S(c). Thus, r(x) can berewritten with such posterior probability distribution functions giventhe independence of x and α and P_(in)=P_(out):${r(x)} = {\frac{{\overset{\_}{P}\left( {{in}❘x} \right)}{\overset{\_}{P}\left( {{in}❘\alpha} \right)}}{{\overset{\_}{P}\left( {{out}❘x} \right)}{\overset{\_}{P}\left( {{out}❘\alpha} \right)}} = \frac{\overset{\_}{P}\left( {{in}❘\left( {x,\alpha} \right)} \right)}{\overset{\_}{P}\left( {{out}❘\left( {x,\alpha} \right)} \right)}}$Thus, this joint likelihood ratio test segmentation is an optimalBayesian binary classification of each voxel with a uniform cost whenthe likelihoods presented herein above are used and the likelihood ratiotest threshold th in is set to one.

FIG. 2 depicts a flow chart of a joint space-intensity likelihood ratiotest based segmentation method, according to an embodiment of theinvention. The segmentation method begins by providing at step 20 animage volume I(x) with a marker x_(p) indicating an approximate locationof a target structure, such as a blob or a tumor. At step 21, a volumeof interest VOI=I(xεΩ), centered at x_(p), is extracted from the imagevolume I. At step 22, an anisotropic Gaussian fitting is performed,resulting in an estimate of the target center u and anisotropic spreadmatrix Σ. At step 23, given the estimated target center and spread (u,Σ) and the VOI=I(xεΩ), the foreground and background spatial andintensity likelihood functions are estimated over a support window S. Atstep 24, for a voxel x in the support window and its associatedintensity α, the likelihood ratio r(x) is computed from thespatial-intensity joint likelihood functions. At step 25, the likelihoodratio test is performed to determine if the voxel is inside or outsidethe target structure: xεin if r(x)≧th, otherwise xεout, where th is athreshold that is optimally set to one, and in and out label the insideand outside of the target structure, respectively. At step 26, thepreceding 2 steps, step 24 and 25, are repeated for all voxels withinthe support window.

A 3D implementation according to an embodiment of the invention wasapplied to delineating a target lung nodule from background lungparenchyma in the presence of other non-target structures such asvessels and lung walls. The performance was evaluated by using highresolution chest CT images of 39 patients including 1310 lung nodules.The images are of size 512×512×400 voxels (depth slightly varies acrossthe patients) with a 12 bit intensity range. For each lung tumor, anapproximate location marker is provided by an expert radiologist. Thesize of VOI is fixed to be 33×33×33 voxels.

FIGS. 3(a)-(d) illustrate examples of segmentation results for fourtumor cases, according to an embodiment of the invention. Each column ofan example corresponds to the segmentation results on yz, xz, xy planes,respectively, passing through the estimated tumor center u. The firstrow of each example depicts the segmentation results without using thederived support window S. In this case, the intensity likelihoods areestimated by using all samples within the 33×3×33 VOI. The second rowdepicts the likelihood ratio segmentation results using the derivedsupport window S. The third row depicts results from a 4Dspace-intensity joint-domain mean shift segmentation. The resultspresented here illustrate that a likelihood ratio based segmentationsolution with a support window successfully performs 3D lung tumorboundary segmentation, while the mean shift and the likelihood ratiowithout S tend to under- and over-estimate the tumor boundary,respectively.

FIGS. 4(a)-(d) shows the intensity likelihood models ƒ(α|in) andƒ(α|out) estimated for the four cases in FIGS. 3(a)-(d), according to anembodiment of the invention. In each example, the dark and light linesdenote the foreground and the background models, respectively. The firstrow of each example illustrate the likelihoods computed with all samplesin the VOI, while the second rows illustrate those computed with thesupport S. The background likelihoods without S, covering a largebackground region, tend to over-sample the lung parenchyma regions withlow intensity values (expressed by a high peak at the low intensityrange). This over-sampling excessively suppresses the backgroundlikelihood at the higher intensity values, resulting in a falsely largeintensity range in which the foreground likelihood surpasses thebackground one. This is an obvious cause for the over-estimation by thissolution. The solution with the support window S, however, employs thebackground intensity model estimated with samples only within S. This ineffect reduces the intensity range to an appropriate size, resulting inbetter segmentation. Example (d) illustrates this effect clearly, wherethe solution without S does not effectively discriminate the intensityinformation.

FIG. 5 illustrates examples of 2D views of 3D segmentation results forfive tumor cases, according to an embodiment of the invention. Thesecross sectional views pass through the estimated tumor center u. Theleft column illustrates the input data. The middle column illustrates ananisotropic Gaussian fitted to the data. A “+” in the image indicatesthe marker x_(p), an “%” indicates the estimated center u, and theellipses in the image indicate the image-plane intersection of the 35%confidence ellipsoid of the estimated Gaussian. The right columnillustrates segmentation results shown as grayscale images with thesegmented regions filled in with a white value. These results illustratethe capability of a joint spatial-intensity likelihood ratio basedsegmentation to handle irregular 3D boundary geometries. The fourth rowof the figure also illustrates a case where the presence of aneighboring lung wall was segmented correctly.

With the 1310 tumor cases, the Gaussian fitting pre-process successfullys approximated the tumor boundary for 1139 cases. Most of the failureswere due to a few isolated voxels near the target boundary being falselysegmented as a part of the target when non-target structures werepresent nearby. This can be mitigated by performing a connectedcomponent analysis as a post-process. After such a post-process, theerror rate reduces to only 1% (11 cases). On average, a method accordingto an embodiment of the invention can run in less than 3 seconds with a2.4 GHz Pentium IV processor, or 3 times faster than a mean shiftsolution.

It is to be understood that the present invention can be implemented invarious forms of hardware, software, firmware, special purposeprocesses, or a combination thereof. In one embodiment, the presentinvention can be implemented in software as an application programtangible embodied on a computer readable program storage device. Theapplication program can be uploaded to, and executed by, a machinecomprising any suitable architecture.

FIG. 6 is a block diagram of an exemplary computer system forimplementing a toboggan-based object characterization scheme accordingto an embodiment of the invention. Referring now to FIG. 6, a computersystem 61 for implementing the present invention can comprise, interalia, a central processing unit (CPU) 62, a memory 63 and aninput/output (I/O) interface 64. The computer system 61 is generallycoupled through the I/O interface 64 to a display 65 and various inputdevices 66 such as a mouse and a keyboard. The support circuits caninclude circuits such as cache, power supplies, clock circuits, and acommunication bus. The memory 63 can include random access memory (RAM),read only memory (ROM), disk drive, tape drive, etc., or a combinationsthereof. The present invention can be implemented as a routine 67 thatis stored in memory 63 and executed by the CPU 62 to process the signalfrom the signal source 68. As such, the computer system 61 is a generalpurpose computer system that becomes a specific purpose computer systemwhen executing the routine 67 of the present invention.

The computer system 61 also includes an operating system and microinstruction code. The various processes and functions described hereincan either be part of the micro instruction code or part of theapplication program (or combination thereof) which is executed via theoperating system. In addition, various other peripheral devices can beconnected to the computer platform such as an additional data storagedevice and a printing device.

It is to be further understood that, because some of the constituentsystem components and method steps depicted in the accompanying figurescan be implemented in software, the actual connections between thesystems components (or the process steps) may differ depending upon themanner in which the present invention is programmed. Given the teachingsof the present invention provided herein, one of ordinary skill in therelated art will be able to contemplate these and similarimplementations or configurations of the present invention.

The particular embodiments disclosed above are illustrative only, as theinvention may be modified and practiced in different but equivalentmanners apparent to those skilled in the art having the benefit of theteachings herein. Furthermore, no limitations are intended to thedetails of construction or design herein shown, other than as describedin the claims below. It is therefore evident that the particularembodiments disclosed above may be altered or modified and all suchvariations are considered within the scope and spirit of the invention.Accordingly, the protection sought herein is as set forth in the claimsbelow.

1. A method for segmenting a digitized image comprising the steps of:providing a digitized volumetric image comprising a plurality ofintensities corresponding to a domain of points in an N-dimensionalspace; providing an approximate location of a target structure in saidimage; estimating a foreground spatial-intensity likelihood functionabout said target structure; estimating a background spatial-intensitylikelihood function about said target structure; and using saidforeground and background spatial-intensity likelihood functions tosegment said target structure by determining whether a point about saidtarget structure is inside said target structure.
 2. The method of claim1, further comprising determining an estimated center and an estimatedspread of said target structure by fitting an N-dimensional anisotropicGaussian function to a volume of interest centered about the approximatelocation and determining the center and the anisotropic spread of saidGaussian function.
 3. The method of claim 2, wherein said foregroundspatial-intensity likelihood function can be factored into a product ofa foreground spatial likelihood function and a foreground intensitylikelihood function, and said background spatial-intensity likelihoodfunction can be factored into a product of a background spatiallikelihood function and a background intensity likelihood function. 4.The method of claim 3, wherein said foreground spatial likelihoodfunction is proportional to said anisotropic Gaussian function, and saidbackground intensity likelihood function is a complement of saidforeground spatial likelihood function.
 5. The method of claim 3,further comprising imposing a window about said target structure,wherein said window is defined as those points whose Mahalanobisdistance from said mean of said Gaussian is less than a predeterminedconstant value, wherein said Mahalanobis distance is computed using saidspread of said Gaussian.
 6. The method of claim 5, wherein said constantvalue is determined by solving${{{{2{\pi\Sigma}}}^{{- 1}/2}{\int_{S{(c)}}{{\exp\left( {{- \frac{1}{2}}\left( {x - u} \right)^{t}{\sum\limits^{- 1}\left( {x - u} \right)}} \right)}{\mathbb{d}x}}}} = \frac{{S(c)}}{2{{2{\pi\Sigma}}}^{1/2}}},$wherein Σ is said spread, c is said constant value, S(c) is said window,x is a point in said window, and u is the center of said targetstructure, and S(c) = ∫_(S(c))𝕕x.
 7. The method of claim 5, whereinsaid foreground intensity likelihood function is proportional to aforeground intensity difference function weighted by said foregroundspatial likelihood function sampled within said window, and saidbackground intensity likelihood function is proportional to a backgroundintensity difference function weighted by said background spatiallikelihood function sampled within said window.
 8. The method of claim7, wherein said proportionality constant is equal to one half the normof the window.
 9. The method of claim 7, wherein the foreground andbackground intensity difference functions comprise Dirac deltafunctions.
 10. The method of claim 7, wherein the foreground andbackground intensity difference functions comprise Parzen functions. 11.The method of claim 1, wherein said step of determining whether a pointabout said target structure is inside said target structure is repeatedfor every point neighboring said target structure to determine whichpoints comprise said target structure.
 12. The method of claim 1,wherein determining whether a point is inside said target structurecomprises comparing a ratio of said foreground and backgroundspatial-intensity likelihood functions calculated at said point to apreset threshold, wherein said point is classified as inside said targetstructure if said ratio is greater than said threshold.
 13. The methodof claim 1, wherein determining whether a point is inside said targetstructure comprises comparing said foreground spatial-intensitylikelihood function ƒ(x, α|in) to said background spatial-intensitylikelihood function ƒ(x, α|out), wherein said point x is classified asinside said target structure if ƒ(x, α|in)>ƒ(x, α|out).
 14. The methodof claim 1, wherein determining whether a point is inside said targetstructure comprises comparing a function F of a ratio of said foregroundlikelihood function to said background likelihood function to F(1) atthe point location, wherein the function F is a member of a family offunctions F:R→R that are monotonically and strictly increasing, whereinsaid point is classified as inside said target structure if the functionof said ratio is greater than F(1).
 15. The method of claim 1, whereindetermining whether a point is inside said target structure comprisescomparing a function F of said foreground likelihood function ƒ(x, α|in)to a function F of said background likelihood function ƒ(x, α|out) atthe point location x, wherein the function F is a member of a family offunctions F:R→R that are monotonically and strictly increasing, andwherein the point x is classified as inside said target structure ifF(ƒ(x, α|in))≧F(ƒ(x, α|out).
 16. A method for segmenting a digitizedimage comprising the steps of: providing a digitized volumetric imagecomprising a plurality of intensities corresponding to a domain ofpoints in an N-dimensional space; identifying a target structure in saidimage; forming a window about said target structure whose size is afunction of the target scale; and performing a jointspace-intensity-likelihood ratio test at each point within said windowto determine whether each said point is within said target structure.17. The method of claim 16, wherein identifying a target structurecomprises providing a marker point within said target structure andfitting an N-dimensional anisotropic Gaussian function about said markerpoint and determining a center and an anisotropic spread of saidGaussian function.
 18. The method of claim 16, wherein performing ajoint space-intensity-likelihood ratio test further comprises:estimating a foreground spatial-intensity likelihood function over awindow encompassing said target structure; estimating a backgroundspatial-intensity likelihood function over said window encompassing saidtarget structure; and comparing a ratio of said foreground andbackground spatial-intensity likelihood functions to a preset thresholdto determine whether a point in said window is inside said targetstructure.
 19. A program storage device readable by a computer, tangiblyembodying a program of instructions executable by the computer toperform the method steps for segmenting a digitized image, the methodcomprising the steps of: providing a digitized volumetric imagecomprising a plurality of intensities corresponding to a domain ofpoints in an N-dimensional space; providing an approximate location of atarget structure in said image; estimating a foregroundspatial-intensity likelihood function about said target structure;estimating a background spatial-intensity likelihood function about saidtarget structure; and using said foreground and backgroundspatial-intensity likelihood functions to segment said target structureby determining whether a point about said target structure is insidesaid target structure.
 20. The computer readable program storage deviceof claim 19, the method further comprising determining an estimatedcenter and an estimated spread of said target structure by fitting anN-dimensional anisotropic Gaussian function to a volume of interestcentered about the approximate location and determining the center andthe anisotropic spread of said Gaussian function.
 21. The computerreadable program storage device of claim 20, wherein said foregroundspatial-intensity likelihood function can be factored into a product ofa foreground spatial likelihood function and a foreground intensitylikelihood function, and said background spatial-intensity likelihoodfunction can be factored into a product of a background spatiallikelihood function and a background intensity likelihood function. 22.The computer readable program storage device of claim 21, wherein saidforeground spatial likelihood function is proportional to saidanisotropic Gaussian function, and said background intensity likelihoodfunction is a complement of said foreground spatial likelihood function.23. The computer readable program storage device of claim 21, the methodfurther comprising imposing a window about said target structure,wherein said window is defined as those points whose Mahalanobisdistance from said mean of said Gaussian is less than a predeterminedconstant value, wherein said Mahalanobis distance is computed using saidspread of said Gaussian.
 24. The computer readable program storagedevice of claim 23, wherein said constant value is determined by solving${{{{2{\pi\Sigma}}}^{{- 1}/2}{\int_{S{(c)}}{{\exp\left( {{- \frac{1}{2}}\left( {x - u} \right)^{t}{\sum\limits^{- 1}\left( {x - u} \right)}} \right)}{\mathbb{d}x}}}} = \frac{{S(c)}}{2{{2{\pi\Sigma}}}^{1/2}}},$wherein Σ is said spread, c is said constant value, S(c) is said window,x is a point in said window, and u is the center of said targetstructure, and S(c) = ∫_(S(c))𝕕x.
 25. The computer readable programstorage device of claim 23, wherein said foreground intensity likelihoodfunction is proportional to a foreground intensity difference functionweighted by said foreground spatial likelihood function sampled withinsaid window, and said background intensity likelihood function isproportional to a background intensity difference function weighted bysaid background spatial likelihood function sampled within said window.26. The computer readable program storage device of claim 25, whereinsaid proportionality constant is equal to one half the norm of thewindow.
 27. The computer readable program storage device of claim 25,wherein the foreground and background intensity difference functionscomprise Dirac delta functions.
 28. The computer readable programstorage device of claim 25, wherein the foreground and backgroundintensity difference functions comprise Parzen functions.
 29. Thecomputer readable program storage device of claim 19, wherein said stepof determining whether a point about said target structure is insidesaid target structure is repeated for every point neighboring saidtarget structure to determine which points comprise said targetstructure.
 30. The computer readable program storage device of claim 19,wherein determining whether a point is inside said target structurecomprises comparing a ratio of said foreground and backgroundspatial-intensity likelihood functions calculated at said point to apreset threshold, wherein said point is classified as inside said targetstructure if said ratio is greater than said threshold.
 31. The computerreadable program storage device of claim 19, wherein determining whethera point is inside said target structure comprises comparing saidforeground spatial-intensity likelihood function ƒ(x, a|in) to saidbackground spatial-intensity likelihood function ƒ(x, α|out), whereinsaid point x is classified as inside said target structure if ƒ(x,α|in)>ƒ(x, α|out).
 32. The computer readable program storage device ofclaim 19, wherein determining whether a point is inside said targetstructure comprises comparing a function F of a ratio of said foregroundlikelihood function to said background likelihood function to F(1) atthe point location, wherein the function F is a member of a family offunctions F:R→R that are monotonically and strictly increasing, whereinsaid point is classified as inside said target structure if the functionof said ratio is greater than F(1).
 33. The computer readable programstorage device of claim 19, wherein determining whether a point isinside said target structure comprises comparing a function F of saidforeground likelihood function ƒ(x, α|in) to a function F of saidbackground likelihood function ƒ(x, α|out) at the point location x,wherein the function F is a member of a family of functions F:R→R thatare monotonically and strictly increasing, and wherein the point x isclassified as inside said target structure if F(ƒ(x, α|in))≧F(ƒ(x,α|out).